//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
//
// --------------------------------------------------------------
//   GEANT 4 - Underground Dark Matter Detector Advanced Example
//
//      For information related to this code contact: Alex Howard
//      e-mail: alexander.howard@cern.ch
// --------------------------------------------------------------
// Comments
//
//                  Underground Advanced
//               by A. Howard and H. Araujo 
//                    (27th November 2001)
//
// DetectorConstruction program
// --------------------------------------------------------------

#include "DMXDetectorConstruction.hh"
#include "DMXDetectorMessenger.hh"

#include "DMXScintSD.hh"
#include "DMXPmtSD.hh"

#include "G4Material.hh"
#include "G4MaterialTable.hh"
#include "G4Element.hh"
#include "G4Isotope.hh"
#include "G4UnitsTable.hh"
#include "G4Box.hh"
#include "G4Tubs.hh"
#include "G4Torus.hh"
#include "G4Sphere.hh"
#include "G4UnionSolid.hh"
#include "G4SubtractionSolid.hh"

#include "G4LogicalVolume.hh"
#include "G4PVPlacement.hh"

#include "MyDeviceParameterisation.hh"
#include "G4PVParameterised.hh"

#include "G4ThreeVector.hh"
#include "G4RotationMatrix.hh"
#include "G4Transform3D.hh"
#include "G4LogicalBorderSurface.hh"
#include "G4LogicalSkinSurface.hh"
#include "G4OpBoundaryProcess.hh"

#include "G4FieldManager.hh"
#include "G4UniformElectricField.hh"
#include "G4TransportationManager.hh"
#include "G4MagIntegratorStepper.hh"
#include "G4EqMagElectricField.hh"
#include "G4ClassicalRK4.hh"
#include "G4ChordFinder.hh"

#include "G4SDManager.hh"

#include "G4VisAttributes.hh"
#include "G4Colour.hh"

#include "G4UserLimits.hh"

#include "G4RunManager.hh"

#include <TROOT.h>
#include <TH1.h>
#include <TH2.h>
#include <TCanvas.h>
#include <TRandom.h>
#include <TNtuple.h>


//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
DMXDetectorConstruction::DMXDetectorConstruction()  
{
  // create commands for interactive definition of time cuts:
  detectorMessenger = new DMXDetectorMessenger(this);

  theUserLimitsForRoom     = 0; 
  theUserLimitsForWorld    = 0; 
  theUserLimitsForDetector = 0; 
  // default time cut = infinite
  //  - note also number of steps cut in stepping action = MaxNoSteps
  theMaxTimeCuts      = DBL_MAX;
  theMaxStepSize      = DBL_MAX;
  theDetectorStepSize = DBL_MAX;
  theWorldStepSize    = 1000.0*mm;
  theRoomTimeCut      = 1000. * nanosecond;
  theMinEkine         = 250.0*eV; // minimum kinetic energy required in volume
  theRoomMinEkine     = 250.0*eV; // minimum kinetic energy required in volume

}


//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
DMXDetectorConstruction::~DMXDetectorConstruction() 
{
  delete theUserLimitsForRoom;
  delete theUserLimitsForWorld;
  delete theUserLimitsForDetector;
  delete detectorMessenger;
}



//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
void DMXDetectorConstruction::DefineMaterials() 
{

#include "DMXDetectorMaterial.icc"

}


/*
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
void DMXDetectorConstruction::DefineField() {

#include "DMXDetectorField.icc"

}
*/


//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
G4VPhysicalVolume* DMXDetectorConstruction::Construct() {

  DefineMaterials();

  // DefineField();

  // make colours
  G4Colour  white   (1.0, 1.0, 1.0) ;
  G4Colour  grey    (0.5, 0.5, 0.5) ;
  G4Colour  lgrey   (.85, .85, .85) ;
  G4Colour  red     (1.0, 0.0, 0.0) ;
  G4Colour  blue    (0.0, 0.0, 1.0) ;
  G4Colour  cyan    (0.0, 1.0, 1.0) ;
  G4Colour  magenta (1.0, 0.0, 1.0) ; 
  G4Colour  purple  (.75, 0.0, 1.0) ; 
  G4Colour  yellow  (1.0, 1.0, 0.0) ;
  G4Colour  orange  (.75, .55, 0.0) ;
  G4Colour  lblue   (0.0, 0.0, .75) ;
  G4Colour  lgreen  (0.0, .75, 0.0) ;
  G4Colour  green   (0.0, 1.0, 0.0) ;
  G4Colour  brown   (0.7, 0.4, 0.1) ;
  //  un-used colours:
  //  G4Colour  black   (0.0, 0.0, 0.0) ;


  // Universe - room wall - CONCRETE ****************************************

  //NB: measured INSIDE of lab, therefore have to add twice wall thickness
  G4double wallThick   = 24.*cm;
  G4double worldWidth  = 470.0*cm + 2.*wallThick; // "x"
  G4double worldLength = 690.0*cm + 2.*wallThick; // "y"
  G4double worldHeight = 280.0*cm + 2.*wallThick; // "z"

  G4Box* world_box = new G4Box("world_box", 0.5*worldWidth, 0.5*worldLength, 0.5*worldHeight );
  world_log  = new G4LogicalVolume(world_box, world_mat, "world_log");
  world_phys = new G4PVPlacement(0, G4ThreeVector(0.,0.,0.), "world_phys", world_log, NULL, false,0);

  //  G4VisAttributes* world_vat= new G4VisAttributes(white);
  world_log->SetVisAttributes(G4VisAttributes::Invisible);
  //world_vat->SetVisibility(true);
  //world_vat->SetVisibility(false);
  //world_log->SetVisAttributes(world_vat);


  // Lab Space - AIR ********************************************************

  G4double labWidth  = worldWidth  - 2.*wallThick; //X
  G4double labLength = worldLength - 2.*wallThick; //Y
  G4double labHeight = worldHeight - 2.*wallThick; //Z

  G4Box* lab_box = new G4Box("lab_box", 0.5*labWidth, 0.5*labLength, 0.5*labHeight );
  lLab  = new G4LogicalVolume(lab_box, lab_mat, "lLab");
  pLab = new G4PVPlacement(0, G4ThreeVector(0.,0.,0.), "pLab", lLab, world_phys, false,0);

  G4VisAttributes* vLab= new G4VisAttributes(white);
  //  lLab->SetVisAttributes(G4VisAttributes::Invisible);
  //  vLab->SetVisibility(true);
  vLab->SetVisibility(false);
  lLab->SetVisAttributes(vLab);

// include room furniture: **************************************************

#include "DMXDetectorRoom.icc"

  //--------Set Visibility
  G4bool IsVisibleOuterJC   = false;
  G4bool IsVisibleInnerJC   = false;
  G4bool IsVisibleGAr       = false;
  G4bool IsVisibleLAr       = false;
  G4bool IsVisibleFid       = false;
  G4bool IsVisibleAnodePad  = false;
  G4bool IsVisibleReflector = false;
  G4bool IsVisibleShaper    = false;
  G4bool IsVisibleGrid      = true;
  G4bool IsVisibleMS        = false;
  G4bool IsVisibleJPMT      = false;
  G4bool IsVisibleMPMT      = false;

  // Now start with detector assembly:
  //*******************J-Chamber**************************// 

  G4double JCHeight = 1050.0*mm;
  
  G4double  FlangeMinR   = 0.;
  G4double  FlangeMaxR   = 170.;
  G4double  TopFlangeHeight = 45.0*mm;
  G4double  ChamSPhi     = 0.*deg;
  G4double  ChamDPhi     = 360.*deg;

  G4double OuterJCHeight = 1095.0*mm;
  G4double OuterJCinnerR = 167.0*mm;
  G4double OuterJCouterR = 170.0*mm;
  G4double InnerJCHeight = 1050.0*mm;
  G4double InnerJCinnerR = 150.0*mm;
  G4double InnerJCouterR = 153.0*mm;

  /*
  G4Tubs* POINT=new G4Tubs("POINT",0.,50.,5.,0.,360.*deg);
  G4LogicalVolume* lPOINT=new G4LogicalVolume(POINT,JC_mat,"lPOINT");
  G4PVPlacement* pPOINT=new G4PVPlacement(0,G4ThreeVector(0,0,-1095.),lPOINT,"pPOINT",lLab,false,0);
  G4VisAttributes* vPOINT = new G4VisAttributes(red);
  vPOINT->SetVisibility(true);
  lPOINT->SetVisAttributes(vPOINT);
  */
  //--- Bottom Flange
  G4double BottomFlangeHeight = 25.0*mm;
  G4double BottomHeight = 45.0*mm;
  G4Tubs* TubsBottomFlange=new G4Tubs("TubsBottomFlange",0.,OuterJCouterR,0.5*BottomFlangeHeight,ChamSPhi,ChamDPhi);
  //--- Outer JCtube
  G4Tubs* TubsOuterJC=new G4Tubs("TubsOuterJC",OuterJCinnerR,OuterJCouterR,0.5*OuterJCHeight,ChamSPhi,ChamDPhi);
  //
  G4UnionSolid* SolidOuterJC=new G4UnionSolid("SolidOuterJC",TubsOuterJC,TubsBottomFlange,G4Transform3D(G4RotationMatrix(),G4ThreeVector(0,0,-0.5*(OuterJCHeight-BottomFlangeHeight))) );

  G4LogicalVolume* lOuterJC=new G4LogicalVolume(SolidOuterJC,JC_mat,"lOuterJC");
  G4PVPlacement* pOuterJC=new G4PVPlacement(0,G4ThreeVector(0,0,-0.5*OuterJCHeight),lOuterJC,"pOuterJC",lLab,false,0);
  G4VisAttributes* vOuterJC = new G4VisAttributes(lgrey);
  vOuterJC->SetVisibility(IsVisibleOuterJC);
  lOuterJC->SetVisAttributes(vOuterJC);


  //--- Inner JCtube
  G4double InnerThick = 2.*mm;
  G4Tubs* TubsInnerJC = new G4Tubs("TubsInnerJC",InnerJCinnerR,InnerJCouterR,0.5*InnerJCHeight,0.*deg,360.*deg);
  G4Tubs* BottomInnerJC = new G4Tubs("BottomInnerJC",0.,InnerJCouterR,0.5*InnerThick*mm,0.*deg,360.*deg);
  G4UnionSolid* SolidInnerJC=new G4UnionSolid("SolidInnerJC",TubsInnerJC,BottomInnerJC,G4Transform3D(G4RotationMatrix(),G4ThreeVector(0,0,-0.5*(InnerJCHeight+InnerThick))) );

  G4LogicalVolume* lInnerJC  = new G4LogicalVolume(SolidInnerJC, JC_mat, "lInnerJC");
  
  G4PVPlacement* pInnerJC = new G4PVPlacement(0, G4ThreeVector(0.,0.,-0.5*InnerJCHeight),lInnerJC,"pInnerJC", lLab, false,0);

  G4VisAttributes* vInnerJC = new G4VisAttributes(lgrey);
  vInnerJC->SetVisibility(IsVisibleInnerJC);
  lInnerJC->SetVisAttributes(vInnerJC);
  
  //top
  //  G4Tubs* solidTop = new G4Tubs("Top",FlangeMinR,FlangeMaxR,0.5*FlangeHeight,ChamSPhi,ChamDPhi);
  //G4LogicalVolume* logicTop = new G4LogicalVolume(solidTop,JC_mat,"Flange");
  // G4PVPlacement* physTop   = new G4PVPlacement(0,G4ThreeVector(0,0,0.5*(FlangeHeight+JCHeight)),logicTop,"physTop",lLab,false,0);
  //G4PVPlacement* physTop   = new G4PVPlacement(0,G4ThreeVector(0,0,0.5*(FlangeHeight+JCHeight)),logicTop,"physTop",lLab,false,0);
  
  //G4VisAttributes* JC_vat = new G4VisAttributes(G4Colour(1.0,1.0,1.0));
  //JC_vat->SetVisibility(true);
  //JC_log->SetVisAttributes(JC_vat);
  
  //--------------------- GAr
  G4double innerGArRadius = 0.*mm;
  G4double outerGArRadius = 150.*mm;
  //  G4double HeightGAr = 355.*mm; //2phase
  G4double HeightGAr = 348.*mm; //1phase
  G4double startAngleGAr    = 0.*deg;
  G4double spanningAngleGAr = 360.*deg;
  G4Tubs* GAr_tube=new G4Tubs("GAr_tube",innerGArRadius,outerGArRadius,0.5*HeightGAr,startAngleGAr,spanningAngleGAr);
  lGAr=new G4LogicalVolume(GAr_tube,GAr_mat,"lGAr"); 
  pGAr=new G4PVPlacement(0,G4ThreeVector(0.,0.,-0.5*HeightGAr),lGAr,"pGAr",lLab,false,0);
  
  G4VisAttributes* vGAr=new G4VisAttributes(cyan);
  vGAr->SetVisibility( IsVisibleGAr );
  lGAr->SetVisAttributes(vGAr);
  
  //--------------------- LAr
  const G4int NUMENTRIES = 3;
  G4double LAr_PP[NUMENTRIES]    = { 9.286*eV, 9.686*eV, 10.186*eV };
  G4double LAr_SCINT[NUMENTRIES] = { 0.170, 1.000, 0.100 };
  G4double LAr_RIND[NUMENTRIES]  = { 1.356, 1.378, 1.418 };
  G4double LAr_ABSL[NUMENTRIES]  = { 1000*m, 1000*m, 1000*m};
  G4double LAr_RAYL[NUMENTRIES]  = { 125.9*cm, 91.07*cm, 56.14*cm}; // Rayleigh scattering length 
  G4MaterialPropertiesTable *LAr_mt = new G4MaterialPropertiesTable();
  LAr_mt->AddProperty("FASTCOMPONENT", LAr_PP, LAr_SCINT, NUMENTRIES);
  LAr_mt->AddProperty("SLOWCOMPONENT", LAr_PP, LAr_SCINT, NUMENTRIES);
  LAr_mt->AddProperty("RINDEX",        LAr_PP, LAr_RIND,  NUMENTRIES);
  LAr_mt->AddProperty("ABSLENGTH",     LAr_PP, LAr_ABSL,  NUMENTRIES);
  LAr_mt->AddProperty("RAYLEIGH",      LAr_PP, LAr_RAYL,  NUMENTRIES);
  LAr_mt->AddConstProperty("SCINTILLATIONYIELD",41667./MeV);
  LAr_mt->AddConstProperty("RESOLUTIONSCALE",1.0);
  LAr_mt->AddConstProperty("FASTTIMECONSTANT",7.*ns);
  LAr_mt->AddConstProperty("SLOWTIMECONSTANT",1600.*ns);
  LAr_mt->AddConstProperty("YIELDRATIO",0.23);

  G4double innerLArRadius = 0.*mm;
  G4double outerLArRadius = 150.*mm;
  G4double HeightLAr = JCHeight - HeightGAr;
  G4double startAngleLAr    = 0.*deg;
  G4double spanningAngleLAr = 360.*deg;
  G4Tubs* LAr_tube=new G4Tubs("LAr_tube",innerLArRadius,outerLArRadius,0.5*HeightLAr,startAngleLAr,spanningAngleLAr);
  lLAr=new G4LogicalVolume(LAr_tube,LAr_mat,"lLAr"); 
  pLAr=new G4PVPlacement(0,G4ThreeVector(0.,0.,-1.*(HeightGAr+0.5*HeightLAr)),lLAr,"pLAr",lLab,false,0);
  G4VisAttributes* vLAr=new G4VisAttributes(blue);
  vLAr->SetVisibility( IsVisibleLAr );
  lLAr->SetVisAttributes(vLAr);

  
  //******************* TPC ********************//
  G4double LArCenterZpos = -1.*(HeightGAr + 0.5*HeightLAr) ;
  G4double GArCenterZpos = -0.5*HeightGAr ;

  //-------------------- Fiducial Volume
  G4double innerFidRadius = 0.*mm;
  G4double outerFidRadius = 128.*mm;
  G4double HeightFid = 190.*mm;
  G4double startAngleFid    = 45.*deg;
  G4double spanningAngleFid = 270.*deg;
  G4double FidBoxSide = outerFidRadius/sqrt(2.);
  G4Tubs* Fid_tube=new G4Tubs("Fid_tube",innerFidRadius,outerFidRadius,0.5*HeightFid,startAngleFid,spanningAngleFid);
  G4Box* Fid_Box=new G4Box("Fid_Box", FidBoxSide, FidBoxSide, 0.5*HeightFid);
  G4UnionSolid* Fid_Sol=new G4UnionSolid("Fid_Sol",Fid_tube,Fid_Box,G4Transform3D(G4RotationMatrix(), G4ThreeVector(0,0,0)));
 
  lFid=new G4LogicalVolume(Fid_Sol,LAr_mat,"lFid"); 
  pFid=new G4PVPlacement(0,G4ThreeVector(0.,0.,-455.*mm-LArCenterZpos),lFid,"pFid",lLAr,false,0);
  G4VisAttributes* vFid=new G4VisAttributes(blue);
  vFid->SetVisibility( IsVisibleFid );
  //  vFid->SetForceSolid(true);
  lFid->SetVisAttributes(vFid);

  //-------- Anode Pad
  G4double APT      =   1.0*mm;
  G4double APRadius = 133.0*mm;  
  G4double APShort  = (133.0/sqrt(2.))*mm; 
  G4double APzpos = -350.*mm;

  G4Tubs* AP_plate1 = new G4Tubs("AP_plate1", 0, APRadius, 0.5*APT, 45.*deg, 270.*deg);
  G4Box* AP_plate2 = new G4Box("AP_plate2", APShort, APShort, 0.5*APT);
  
  G4UnionSolid* AP_sol = new G4UnionSolid("AP_sol", AP_plate1, AP_plate2, G4Transform3D(G4RotationMatrix(), G4ThreeVector(0,0,0)));

  lAnode = new G4LogicalVolume(AP_sol, Anode_mat, "AP_log");
  pAnode = new G4PVPlacement(0,G4ThreeVector(0., 0., APzpos-LArCenterZpos),lAnode,"lAnode",lLAr,false,0);//in GAr
  //  pAnode = new G4PVPlacement(0,G4ThreeVector(0., 0., APzpos-GArCenterZpos),lAnode,"lAnode",lGAr,false,0);//in LAr
  
  G4VisAttributes* AP_vat = new G4VisAttributes(orange);
  AP_vat->SetVisibility(IsVisibleAnodePad);
  //  AP_vat->SetForceSolid(true);
  lAnode->SetVisAttributes(AP_vat);
  
  //-------------- Reflector WLS
  /*
  G4double RefWLSRout   = 128.*mm;
  G4double RefWLSRin    = 127.8*mm;
  G4double RefWLSHeight = 180.*mm;
  G4double RefWLSZpos   = -455.*mm; //Mid of AG and CG

  G4Tubs* RefWLSTubs=new G4Tubs("RefWLSTubs",RefWLSRin,RefWLSRout,0.5*RefWLSHeight,45.*deg,270.*deg);
  G4Box* RefWLSBox=new G4Box("RefWLSBox",(RefWLSRout-RefWLSRin)*cos(45.*deg),RefWLSRout*sin(45.*deg),0.5*RefWLSHeight );
  G4UnionSolid* MRefWLS=new G4UnionSolid("MRefWLS",RefWLSTubs,RefWLSBox,0,G4ThreeVector(0.5*(RefWLSRout+RefWLSRin)*cos(45.*deg),0.,0.));

  G4LogicalVolume* LRefWLS=new G4LogicalVolume(MRefWLS,WLS_mat,"LRefWLS");
  G4PVPlacement* PRefWLS=new G4PVPlacement(0,G4ThreeVector(0.,0.,RefWLSZpos-LArCenterZpos),LRefWLS,"PRefWLS",lLAr,false,0);
  G4VisAttributes* VRefWLS=new G4VisAttributes(blue);
  VRefWLS->SetVisibility( IsVisibleReflector );
  LRefWLS->SetVisAttributes( VRefWLS );
  */
  //-------------- Reflector
  /*
  G4double ReflectorRout   = 129.*mm;
  G4double ReflectorRin    = 128.*mm;
  G4double ReflectorHeight = 180.*mm;
  G4double ReflectorZpos   = -455.*mm; //Mid of AG and CG

  G4Tubs* RefTubs=new G4Tubs("RefTubs",ReflectorRin,ReflectorRout,0.5*ReflectorHeight,45.*deg,270.*deg);
  G4Box* RefBox=new G4Box("RefBox",(ReflectorRout-ReflectorRin)*cos(45.*deg),ReflectorRout*sin(45.*deg),0.5*ReflectorHeight );
  G4UnionSolid* MReflector=new G4UnionSolid("MReflector",RefTubs,RefBox,0,G4ThreeVector(0.5*(ReflectorRout+ReflectorRin)*cos(45.*deg),0.,0.));

  G4LogicalVolume* LReflector=new G4LogicalVolume(MReflector,PTFE_mat,"LReflector");
  G4PVPlacement* PReflector=new G4PVPlacement(0,G4ThreeVector(0.,0.,ReflectorZpos-LArCenterZpos),LReflector,"PReflector",lLAr,false,0);
  G4VisAttributes* VReflector=new G4VisAttributes(grey);
  VReflector->SetVisibility( IsVisibleReflector );
  LReflector->SetVisAttributes( VReflector );
  */

  //-------------- Shaper
  G4RotationMatrix* RotM = new G4RotationMatrix;
  RotM->rotateX(90.*deg);

  const G4int nShaper = 18;
  G4double ShaperRmin = 0.5*3. *mm;
  G4double ShaperRmax = 0.5*3.2  *mm;
  G4double ShaperRtor = (133.-ShaperRmax)*mm;
  G4double ShaperSphi = 45. *deg;
  G4double ShaperDphi = 270.*deg;
  G4double ShaperFirstZpos = -370.*mm;

  G4Torus* MShaper1=new G4Torus("MShaper1",ShaperRmin,ShaperRmax,ShaperRtor,ShaperSphi,ShaperDphi);
  G4Tubs* MShaper2=new G4Tubs("MShaper2",ShaperRmin,ShaperRmax, 0.5*(ShaperRtor+ShaperRmax)*sqrt(2), 0.*deg, 360.*deg);
  G4UnionSolid* MShaper=new G4UnionSolid("MShaper",MShaper1,MShaper2,RotM,G4ThreeVector(ShaperRtor/sqrt(2),0,0));

  G4LogicalVolume* LShaper[nShaper];
  G4PVPlacement*   PShaper[nShaper];
  G4VisAttributes* VShaper=new G4VisAttributes(lgrey);;
  VShaper->SetVisibility( IsVisibleShaper );

  for(int ishaper=0;ishaper<nShaper;ishaper++){
    G4String lname = Form("LShaper_%d",ishaper);
    G4String pname = Form("PShaper_%d",ishaper);
    G4ThreeVector posshaper(0.*mm,0.*mm, ShaperFirstZpos - LArCenterZpos -10.*ishaper*mm);
    LShaper[ishaper]=new G4LogicalVolume(MShaper,JC_mat,lname);
    PShaper[ishaper]=new G4PVPlacement(0,posshaper,LShaper[ishaper],pname,lLAr,false,0);
    LShaper[ishaper]->SetVisAttributes( VShaper );
  }

  //-------------------------------Grid
  G4RotationMatrix* RotMg = new G4RotationMatrix;
  RotMg->rotateZ(90.*deg);
  G4double GridThick = 3.1*mm;
  //  G4double GridThick = 1.55*mm;
  G4double GridRout = 133.*mm;
  G4double GridRin  = 123.*mm;
  G4double WireWidth  = 0.1*mm;
  G4double WireSpace  = 4.0*mm;

  G4Tubs* tGrid1=new G4Tubs("tGrid1",GridRin,GridRout, 0.5*GridThick, (45.-1.81)*deg, (270.+1.81)*deg);
  G4Box* tGrid2=new G4Box("tGrid2",0.5*182.04*mm,0.5*10.*mm,0.5*GridThick);
  G4UnionSolid* grid=new G4UnionSolid("grid",tGrid1,tGrid2,RotMg,G4ThreeVector(91.97*mm,0,0));  

  const G4int ngrid=4;
  const G4int nwire=31;
  G4double grid_z[ngrid]={(-360.-LArCenterZpos)*mm,
			  (-550.-LArCenterZpos)*mm,
			  (-570.-LArCenterZpos)*mm,
			  (-590.-LArCenterZpos)*mm};

  G4Box* twire1[ngrid][nwire];
  G4LogicalVolume* Lwire1[ngrid][nwire];
  G4PVPlacement*   Pwire1[ngrid][nwire];
  G4Box* twire2[ngrid][nwire];
  G4LogicalVolume* Lwire2[ngrid][nwire];
  G4PVPlacement*   Pwire2[ngrid][nwire];
  G4Box* twire3[ngrid][nwire];
  G4LogicalVolume* Lwire3[ngrid][nwire];
  G4PVPlacement*   Pwire3[ngrid][nwire];
  G4Box* twire4[ngrid][nwire];
  G4LogicalVolume* Lwire4[ngrid][nwire];
  G4PVPlacement*   Pwire4[ngrid][nwire];

  G4LogicalVolume* LGrid[ngrid];
  G4PVPlacement*   PGrid[ngrid];
  G4VisAttributes* VGrid=new G4VisAttributes(white);
  VGrid->SetVisibility( IsVisibleGrid );
  
  for(int igrid=0;igrid<ngrid;igrid++){
    LGrid[igrid]=new G4LogicalVolume(grid,JC_mat,Form("LGrid%d",igrid));
    LGrid[igrid]->SetVisAttributes( VGrid );
    PGrid[igrid]=new G4PVPlacement(0,G4ThreeVector(0,0,grid_z[igrid]),LGrid[igrid],
				   Form("PGrid%d",igrid),lLAr,false,0);
    //1
    for(int iwire=1;iwire<22;iwire++){
      G4double length=sqrt(GridRin*GridRin-pow(WireSpace*iwire,2));
      G4String wname = Form("iwire%d_%d",igrid,iwire);
      G4String lname = Form("Lwire%d_%d",igrid,iwire);
      G4String pname = Form("Pwire%d_%d",igrid,iwire);
      twire1[igrid][iwire]=new G4Box(wname,WireWidth*0.5,length,WireWidth*0.5);
      Lwire1[igrid][iwire]=new G4LogicalVolume(twire1[igrid][iwire],JC_mat,lname);
      Lwire1[igrid][iwire]->SetVisAttributes( VGrid );
      Pwire1[igrid][iwire]=new G4PVPlacement(0,G4ThreeVector(WireSpace*iwire,0,grid_z[igrid]),Lwire1[igrid][iwire],pname,lLAr,false,0);
    }
    //2
    for(int iwire=0;iwire<nwire;iwire++){
      G4double length=sqrt(GridRin*GridRin-pow(WireSpace*iwire,2));
      G4String wname = Form("iwire%d_%d",igrid,iwire);
      G4String lname = Form("Lwire%d_%d",igrid,iwire);
      G4String pname = Form("Pwire%d_%d",igrid,iwire);
      twire2[igrid][iwire]=new G4Box(wname,WireWidth*0.5,length,WireWidth*0.5);
      Lwire2[igrid][iwire]=new G4LogicalVolume(twire2[igrid][iwire],JC_mat,lname);
      Lwire2[igrid][iwire]->SetVisAttributes( VGrid );
      Pwire2[igrid][iwire]=new G4PVPlacement(0,G4ThreeVector(-1*WireSpace*iwire,0,grid_z[igrid]),Lwire2[igrid][iwire],pname,lLAr,false,0);
    }
    //3
    for(int iwire=1;iwire<nwire;iwire++){
      G4double xpos = 0.;
      G4double length=sqrt(GridRin*GridRin-pow(WireSpace*iwire,2));
      if(iwire<22){
	xpos = 0.5*(-length+86.97);
	length=(sqrt(GridRin*GridRin-pow(WireSpace*iwire,2))+86.97*mm)*0.5;
      }
      G4String wname = Form("iwire%d_%d",igrid,iwire);
      G4String lname = Form("Lwire%d_%d",igrid,iwire);
      G4String pname = Form("Pwire%d_%d",igrid,iwire);
      twire3[igrid][iwire]=new G4Box(wname,length,WireWidth*0.5,WireWidth*0.5);
      Lwire3[igrid][iwire]=new G4LogicalVolume(twire3[igrid][iwire],JC_mat,lname);
      Lwire3[igrid][iwire]->SetVisAttributes( VGrid );
      Pwire3[igrid][iwire]=new G4PVPlacement(0,
				      G4ThreeVector(xpos*mm,
						    WireSpace*iwire,
						    grid_z[igrid]),
				      Lwire3[igrid][iwire],
				      pname,lLAr,false,0);
    }
    //4
    for(int iwire=0;iwire<nwire;iwire++){
      G4double xpos = 0.;
      G4double length=sqrt(GridRin*GridRin-pow(WireSpace*iwire,2));
      if(iwire<22){
	xpos = 0.5*(-length+86.97);
	length=(sqrt(GridRin*GridRin-pow(WireSpace*iwire,2))+86.97*mm)*0.5;
      }
      G4String wname = Form("iwire%d_%d",igrid,iwire);
      G4String lname = Form("Lwire%d_%d",igrid,iwire);
      G4String pname = Form("Pwire%d_%d",igrid,iwire);
      twire4[igrid][iwire]=new G4Box(wname,length,WireWidth*0.5,WireWidth*0.5);
      Lwire4[igrid][iwire]=new G4LogicalVolume(twire4[igrid][iwire],JC_mat,lname);
      Lwire4[igrid][iwire]->SetVisAttributes( VGrid );
      Pwire4[igrid][iwire]=new G4PVPlacement(0,
				      G4ThreeVector(xpos*mm,
						    -1*WireSpace*iwire,
						    grid_z[igrid]),
				      Lwire4[igrid][iwire],
				      pname,lLAr,false,0);
    }
  }

  //---------- Metal Sheet (for purity)

  G4double MSHeight =   1.0*mm;
  G4double MSSide  = 30.*mm; 
  //  G4double MSxpos = 0.*mm;
  G4double MSypos = -80.*mm;
  //  G4double MSypos = -80.*mm;//0.*mm;
  G4double MSxpos = 0.*mm;
  G4double MSzpos = -569.5*mm;

  G4Box* MS_box = new G4Box("MS_box", 0.5*MSSide, 0.5*MSSide, 0.5*MSHeight);
  
  lMS = new G4LogicalVolume(MS_box, Anode_mat, "MS_log");
  pMS = new G4PVPlacement(0,G4ThreeVector(MSxpos,MSypos,MSzpos-LArCenterZpos),lMS,"lMS",lLAr,false,0);
  
  G4VisAttributes* vMS = new G4VisAttributes(orange);
  vMS->SetVisibility(IsVisibleMS);
  lMS->SetVisAttributes(vMS);

  //---------------------------------H1161 PMT
  //_____ J-PMT Window
  /*
  G4double JPMT_WRadius = 38.0*mm;
  G4double JPMT_WHeight = 2.*mm; 
  G4ThreeVector JPMT_Wpos(0.,0.,-601.*mm - LArCenterZpos);

  G4Tubs* JPMT_Wtubs = new G4Tubs("JPMT_Wtubs",0.,JPMT_WRadius,0.5*JPMT_WHeight,0.*deg,360.*deg);

  lJPMT_W = new G4LogicalVolume(JPMT_Wtubs, glass_mat, "lJPMT_W");
  pJPMT_W = new G4PVPlacement(0,JPMT_Wpos,"pJPMT_W", lJPMT_W, pLAr, false,0);

  G4VisAttributes* vJPMT_W = new G4VisAttributes(purple);
  vJPMT_W->SetVisibility( IsVisibleJPMT );
  lJPMT_W->SetVisAttributes(vJPMT_W);

  //_____ WLS for J-PMT
  G4double JPMT_WLSRadius = 38.0*mm;
  G4double JPMT_WLSHeight = 0.2*mm; 
  G4ThreeVector JPMT_WLSpos(0.,0.,-599.9*mm - LArCenterZpos);

  G4Tubs* JPMT_WLStubs = new G4Tubs("JPMT_WLStubs",0., JPMT_WLSRadius, 0.5*JPMT_WLSHeight, 0.*deg, 360.*deg);

  lJPMT_WLS = new G4LogicalVolume(JPMT_WLStubs, WLS_mat, "lJPMT_WLS");
  pJPMT_WLS = new G4PVPlacement(0,JPMT_WLSpos,"pJPMT_WLS", lJPMT_WLS, pLAr, false,0);

  G4VisAttributes* vJPMT_WLS = new G4VisAttributes(blue);
  vJPMT_WLS->SetVisibility( IsVisibleJPMT );
  lJPMT_WLS->SetVisAttributes(vJPMT_WLS);
  */
  //_____ J-PMT tube
  G4double JPMT_TUBERadiusOuter = 38.75*mm;
  G4double JPMT_TUBERadiusInner = 38.0*mm;
  G4double JPMT_TUBEHeight = 123.*mm;
  G4ThreeVector JPMT_TUBEpos(0.,0.,-600.*mm-0.5*JPMT_TUBEHeight - LArCenterZpos);
  G4Tubs* JPMT_TUBEtubs = new G4Tubs("JPMT_TUBEtubs",JPMT_TUBERadiusInner, JPMT_TUBERadiusOuter, 0.5*JPMT_TUBEHeight, 0.*deg, 360.*deg);

  lJPMT_TUBE = new G4LogicalVolume(JPMT_TUBEtubs, PTFE_mat, "lJPMT_TUBE");
  pJPMT_TUBE = new G4PVPlacement(0,JPMT_TUBEpos,"pJPMT_TUBE", lJPMT_TUBE, pLAr, false,0);

  G4VisAttributes* vJPMT_TUBE = new G4VisAttributes(white);
  //vJPMT_TUBE->SetForceSolid(true);
  vJPMT_TUBE->SetVisibility( IsVisibleJPMT );
  lJPMT_TUBE->SetVisAttributes(vJPMT_TUBE);

  //_____ Photo-Cahode
  G4double JPMT_PCRadius = 32.0*mm;
  G4double JPMT_PCHeight = 0.2*mm; 
  G4ThreeVector JPMT_PCpos(0.,0.,-603.*mm - LArCenterZpos);

  G4Tubs* JPMT_PCtubs = new G4Tubs("JPMT_PCtubs",0., JPMT_PCRadius, 0.5*JPMT_PCHeight, 0.*deg, 360.*deg);

  lJPMT_PC = new G4LogicalVolume(JPMT_PCtubs, Bialkali_mat, "lJPMT_PC");
  pJPMT_PC = new G4PVPlacement(0,JPMT_PCpos,"pJPMT_PC", lJPMT_PC, pLAr, false,0);

  G4VisAttributes* vJPMT_PC = new G4VisAttributes(brown);
  //JPMT_vat->SetForceSolid(true);
  vJPMT_PC->SetVisibility( IsVisibleJPMT );
  lJPMT_PC->SetVisAttributes(vJPMT_PC);

  //-----------------------------------MEG PMT Photo-Cahode
  G4double MPMTFromCenter = 80.*mm;
  G4double MPMT1ang =  45.*deg;
  G4double MPMT2ang = 135.*deg;
  G4double MPMT3ang = 225.*deg;
  G4double MPMT4ang = 315.*deg;

  //_____ Meg-PMT Window
  /*
  G4double MPMT_WRadius = 27.8*mm;
  G4double MPMT_WHeight = 2.*mm; 
  G4double MPMT_Wzpos = -601.*mm;
  G4Tubs* MPMT_Wtubs = new G4Tubs("MPMT_Wtubs",0.,MPMT_WRadius,0.5*MPMT_WHeight,0.*deg,360.*deg);

  lMPMT_W = new G4LogicalVolume(MPMT_Wtubs, glass_mat, "lMPMT_W");

  G4ThreeVector MPMT1_Wpos(MPMTFromCenter*cos(MPMT1ang),MPMTFromCenter*sin(MPMT1ang),MPMT_Wzpos-LArCenterZpos);
  G4ThreeVector MPMT2_Wpos(MPMTFromCenter*cos(MPMT2ang),MPMTFromCenter*sin(MPMT2ang),MPMT_Wzpos-LArCenterZpos);
  G4ThreeVector MPMT3_Wpos(MPMTFromCenter*cos(MPMT3ang),MPMTFromCenter*sin(MPMT3ang),MPMT_Wzpos-LArCenterZpos);
  G4ThreeVector MPMT4_Wpos(MPMTFromCenter*cos(MPMT4ang),MPMTFromCenter*sin(MPMT4ang),MPMT_Wzpos-LArCenterZpos);
  pMPMT1_W = new G4PVPlacement(0,MPMT1_Wpos,"pMPMT1_W", lMPMT_W, pLAr, false,0);
  pMPMT2_W = new G4PVPlacement(0,MPMT2_Wpos,"pMPMT2_W", lMPMT_W, pLAr, false,0);
  pMPMT3_W = new G4PVPlacement(0,MPMT3_Wpos,"pMPMT3_W", lMPMT_W, pLAr, false,0);
  pMPMT4_W = new G4PVPlacement(0,MPMT4_Wpos,"pMPMT4_W", lMPMT_W, pLAr, false,0);

  G4VisAttributes* vMPMT_W = new G4VisAttributes(magenta);
  vMPMT_W->SetVisibility( IsVisibleMPMT );
  lMPMT_W->SetVisAttributes(vMPMT_W);

  //_____ Meg-PMT WLS
  G4double MPMT_WLSRadius = 27.8*mm;
  G4double MPMT_WLSHeight = 0.2*mm;
  G4double MPMT_WLSzpos = -599.9*mm; 
  G4Tubs* MPMT_WLStubs = new G4Tubs("MPMT_WLStubs",0.,MPMT_WLSRadius,0.5*MPMT_WLSHeight,0.*deg,360.*deg);

  lMPMT_WLS = new G4LogicalVolume(MPMT_WLStubs, WLS_mat, "lMPMT_WLS");

  G4ThreeVector MPMT1_WLSpos(MPMTFromCenter*cos(MPMT1ang),MPMTFromCenter*sin(MPMT1ang),MPMT_WLSzpos-LArCenterZpos);
  G4ThreeVector MPMT2_WLSpos(MPMTFromCenter*cos(MPMT2ang),MPMTFromCenter*sin(MPMT2ang),MPMT_WLSzpos-LArCenterZpos);
  G4ThreeVector MPMT3_WLSpos(MPMTFromCenter*cos(MPMT3ang),MPMTFromCenter*sin(MPMT3ang),MPMT_WLSzpos-LArCenterZpos);
  G4ThreeVector MPMT4_WLSpos(MPMTFromCenter*cos(MPMT4ang),MPMTFromCenter*sin(MPMT4ang),MPMT_WLSzpos-LArCenterZpos);
  pMPMT1_WLS = new G4PVPlacement(0,MPMT1_WLSpos,"pMPMT1_WLS", lMPMT_WLS, pLAr, false,0);
  pMPMT2_WLS = new G4PVPlacement(0,MPMT2_WLSpos,"pMPMT2_WLS", lMPMT_WLS, pLAr, false,0);
  pMPMT3_WLS = new G4PVPlacement(0,MPMT3_WLSpos,"pMPMT3_WLS", lMPMT_WLS, pLAr, false,0);
  pMPMT4_WLS = new G4PVPlacement(0,MPMT4_WLSpos,"pMPMT4_WLS", lMPMT_WLS, pLAr, false,0);

  G4VisAttributes* vMPMT_WLS = new G4VisAttributes(blue);
  vMPMT_WLS->SetVisibility( IsVisibleMPMT );
  lMPMT_WLS->SetVisAttributes(vMPMT_WLS);
  */
  //_____ Meg-PMT Tube
  G4double MPMT_TUBERadiusOuter = 28.5*mm;
  G4double MPMT_TUBERadiusInner = 27.8*mm;
  G4double MPMT_TUBEHeight = 32.5*mm;
  G4double MPMT_TUBEzpos = -600.*mm-0.5*MPMT_TUBEHeight; 
  G4Tubs* MPMT_TUBEtubs = new G4Tubs("MPMT_TUBEtubs",MPMT_TUBERadiusInner,MPMT_TUBERadiusOuter,0.5*MPMT_TUBEHeight,0.*deg,360.*deg);

  lMPMT_TUBE = new G4LogicalVolume(MPMT_TUBEtubs, PTFE_mat, "lMPMT_WLS");

  G4ThreeVector MPMT1_TUBEpos(MPMTFromCenter*cos(MPMT1ang),MPMTFromCenter*sin(MPMT1ang),MPMT_TUBEzpos-LArCenterZpos);
  G4ThreeVector MPMT2_TUBEpos(MPMTFromCenter*cos(MPMT2ang),MPMTFromCenter*sin(MPMT2ang),MPMT_TUBEzpos-LArCenterZpos);
  G4ThreeVector MPMT3_TUBEpos(MPMTFromCenter*cos(MPMT3ang),MPMTFromCenter*sin(MPMT3ang),MPMT_TUBEzpos-LArCenterZpos);
  G4ThreeVector MPMT4_TUBEpos(MPMTFromCenter*cos(MPMT4ang),MPMTFromCenter*sin(MPMT4ang),MPMT_TUBEzpos-LArCenterZpos);
  pMPMT1_TUBE = new G4PVPlacement(0,MPMT1_TUBEpos,"pMPMT1_TUBE", lMPMT_TUBE, pLAr, false,0);
  pMPMT2_TUBE = new G4PVPlacement(0,MPMT2_TUBEpos,"pMPMT2_TUBE", lMPMT_TUBE, pLAr, false,0);
  pMPMT3_TUBE = new G4PVPlacement(0,MPMT3_TUBEpos,"pMPMT3_TUBE", lMPMT_TUBE, pLAr, false,0);
  pMPMT4_TUBE = new G4PVPlacement(0,MPMT4_TUBEpos,"pMPMT4_TUBE", lMPMT_TUBE, pLAr, false,0);

  G4VisAttributes* vMPMT_TUBE = new G4VisAttributes(white);
  vMPMT_TUBE->SetVisibility( IsVisibleMPMT );
  lMPMT_TUBE->SetVisAttributes(vMPMT_TUBE);

  //_____ Meg-PMT Photo-Cathode
  G4double MPMT_PCRadius = 23.0*mm;
  G4double MPMT_PCHeight = 0.2*mm;
  G4double MPMT_PCzpos = -603.*mm; 
  G4Tubs* MPMT_PCtubs = new G4Tubs("MPMT_PCtubs",0.*cm, MPMT_PCRadius, 0.5*MPMT_PCHeight, 0.*deg, 360.*deg);
  G4VisAttributes* vMPMT_PC = new G4VisAttributes(brown);
  //MPMT_vat->SetForceSolid(true);
  vMPMT_PC->SetVisibility( IsVisibleMPMT );

  //+++MPMT1
  G4ThreeVector MPMT1_PCpos(MPMTFromCenter*cos(MPMT1ang),MPMTFromCenter*sin(MPMT1ang),MPMT_PCzpos-LArCenterZpos);
  lMPMT1_PC = new G4LogicalVolume(MPMT_PCtubs, Bialkali_mat, "lMPMT1_PC");
  pMPMT1_PC = new G4PVPlacement(0,MPMT1_PCpos,"pMPMT1_PC", lMPMT1_PC, pLAr, false,0);
  lMPMT1_PC->SetVisAttributes(vMPMT_PC);
  //+++MPMT2
  G4ThreeVector MPMT2_PCpos(MPMTFromCenter*cos(MPMT2ang),MPMTFromCenter*sin(MPMT2ang),MPMT_PCzpos-LArCenterZpos);
  lMPMT2_PC = new G4LogicalVolume(MPMT_PCtubs, Bialkali_mat, "lMPMT2_PC");
  pMPMT2_PC = new G4PVPlacement(0,MPMT2_PCpos,"pMPMT2_PC", lMPMT2_PC, pLAr, false,0);
  lMPMT2_PC->SetVisAttributes(vMPMT_PC);
  //+++MPMT3
  G4ThreeVector MPMT3_PCpos(MPMTFromCenter*cos(MPMT3ang),MPMTFromCenter*sin(MPMT3ang),MPMT_PCzpos-LArCenterZpos);
  lMPMT3_PC = new G4LogicalVolume(MPMT_PCtubs, Bialkali_mat, "lMPMT3_PC");
  pMPMT3_PC = new G4PVPlacement(0,MPMT3_PCpos,"pMPMT3_PC", lMPMT3_PC, pLAr, false,0);
  lMPMT3_PC->SetVisAttributes(vMPMT_PC);
  //+++MPMT4
  G4ThreeVector MPMT4_PCpos(MPMTFromCenter*cos(MPMT4ang),MPMTFromCenter*sin(MPMT4ang),MPMT_PCzpos-LArCenterZpos);
  lMPMT4_PC = new G4LogicalVolume(MPMT_PCtubs, Bialkali_mat, "lMPMT4_PC");
  pMPMT4_PC = new G4PVPlacement(0,MPMT4_PCpos,"pMPMT4_PC", lMPMT4_PC, pLAr, false,0);
  lMPMT4_PC->SetVisAttributes(vMPMT_PC);


  //********************** Surface
  static const G4double LambdaE = 2.0 * 3.14159265358979323846 * 1.973269602e-16*m*GeV;
  G4double LowOPWaveLength  = 115.0*nanometer;
  G4double HighOPWaveLength = 598.0*nanometer;
  G4double PPCKOVHighE = LambdaE / LowOPWaveLength;
  G4double PPCKOVLowE = LambdaE / HighOPWaveLength;
  //G4cout<<"+++++++ LowE = "<<PPCKOVLowE / eV <<G4endl;
  //G4cout<<"+++++++HighE = "<<PPCKOVHighE/ eV <<G4endl;

  //-------- JChamber-LAr surface
  G4OpticalSurface * OpJCSurface = new G4OpticalSurface
    ("JCSurface", unified, ground, dielectric_metal);
  G4LogicalBorderSurface* JCSurface;
  JCSurface = new G4LogicalBorderSurface("JChamber", pLAr, pInnerJC, OpJCSurface);

  G4double JC_PP[2] = { PPCKOVLowE,PPCKOVHighE};
  G4double JC_REFL[2] = { 0.0, 0.0 };
  //  G4double JC_REFL[2] = { 0.1, 0.1 };
  G4MaterialPropertiesTable *JC_mt = new G4MaterialPropertiesTable();
  JC_mt->AddProperty("REFLECTIVITY", JC_PP, JC_REFL, 2);
  OpJCSurface->SetMaterialPropertiesTable(JC_mt);

  //-------- GAr-LAr Surface
  G4OpticalSurface * OpGArLArSurface = new G4OpticalSurface
    ("GArLArSurface", glisur, polished, dielectric_dielectric);
  G4LogicalBorderSurface* GArLArSurface;
  GArLArSurface = new G4LogicalBorderSurface("GArLArSurface", pLAr, pGAr, OpGArLArSurface);
  OpGArLArSurface->SetMaterialPropertiesTable(LAr_mt);

  //-------- Anode Pad surface
  G4OpticalSurface * OpAPSurface = new G4OpticalSurface
    ("APSurface", unified, ground, dielectric_metal);
  G4LogicalBorderSurface* APSurface;
  APSurface = new G4LogicalBorderSurface("Anode",pLAr,pAnode, OpAPSurface);

  G4double AP_PP[3] = { 2.00*eV, 9.00*eV, 10.00*eV };
  G4double AP_REFL[3] = { 0.0, 0.0, 0.0 };
  //  G4double AP_REFL[3] = { 0.1, 0.1, 0.1 };
  G4MaterialPropertiesTable *AP_mt = new G4MaterialPropertiesTable();
  AP_mt->AddProperty("REFLECTIVITY", AP_PP, AP_REFL, 3);
  OpAPSurface->SetMaterialPropertiesTable(AP_mt);

  //-------- Reflector surface
  /*
  G4OpticalSurface * OpReflectorSurface = new G4OpticalSurface
    ("ReflectorSurface", unified, ground, dielectric_dielectric);
  G4LogicalBorderSurface* ReflectorSurface=new G4LogicalBorderSurface("Reflector",pLAr,PReflector, OpReflectorSurface);
  //G4LogicalSkinSurface* ReflectorSurface = new G4LogicalSkinSurface("Reflector",LReflector, OpReflectorSurface);

  G4double Reflector_PP[3] = { 2.00*eV, 9.00*eV, 10.00*eV };
  G4double Reflector_REFL[3] = { 0.9, 0.9, 0.9 };
  G4MaterialPropertiesTable *Reflector_mt = new G4MaterialPropertiesTable();
  Reflector_mt->AddProperty("REFLECTIVITY", Reflector_PP, Reflector_REFL, 3);
  OpReflectorSurface->SetMaterialPropertiesTable(Reflector_mt);
  */
  //------- Shaper surface
  G4OpticalSurface * OpShaperSurface[nShaper];
  G4LogicalBorderSurface* ShaperSurface;
  for(int ishaper=0;ishaper<nShaper;ishaper++){
    OpShaperSurface[ishaper] = new G4OpticalSurface
      ("ShaperSurface", unified, ground, dielectric_metal);
    ShaperSurface = new G4LogicalBorderSurface(Form("ShaperSur_%d",ishaper), pLAr, PShaper[ishaper], OpShaperSurface[ishaper]);
    OpShaperSurface[ishaper]->SetMaterialPropertiesTable(JC_mt);
  }


  //------- Grid Surface
  G4OpticalSurface * OpGridSurface[ngrid];
  G4OpticalSurface * OpWireSurface1[ngrid][nwire];
  G4OpticalSurface * OpWireSurface2[ngrid][nwire];
  G4OpticalSurface * OpWireSurface3[ngrid][nwire];
  G4OpticalSurface * OpWireSurface4[ngrid][nwire];
  G4LogicalBorderSurface* GridSurface;
  G4LogicalBorderSurface* WireSurface1;
  G4LogicalBorderSurface* WireSurface2;
  G4LogicalBorderSurface* WireSurface3;
  G4LogicalBorderSurface* WireSurface4;
  for(int igrid=0;igrid<ngrid;igrid++){
    OpGridSurface[igrid] = new G4OpticalSurface(Form("GridSurface_%d",igrid), unified, ground, dielectric_metal);
    GridSurface = new G4LogicalBorderSurface("GridSur", pLAr, PGrid[igrid], OpGridSurface[igrid]);
    OpGridSurface[igrid]->SetMaterialPropertiesTable(JC_mt);
    for(int iwire=0;iwire<nwire;iwire++){
      OpWireSurface1[igrid][iwire] = new G4OpticalSurface(Form("WireSurface1_%d_%d",igrid,iwire), unified, ground, dielectric_metal);
      OpWireSurface2[igrid][iwire] = new G4OpticalSurface(Form("WireSurface2_%d_%d",igrid,iwire), unified, ground, dielectric_metal);
      OpWireSurface3[igrid][iwire] = new G4OpticalSurface(Form("WireSurface3_%d_%d",igrid,iwire), unified, ground, dielectric_metal);
      OpWireSurface4[igrid][iwire] = new G4OpticalSurface(Form("WireSurface4_%d_%d",igrid,iwire), unified, ground, dielectric_metal);
      WireSurface1 = new G4LogicalBorderSurface("WireSur1",pLAr,Pwire1[igrid][iwire],OpWireSurface1[igrid][iwire]); 
      WireSurface2 = new G4LogicalBorderSurface("WireSur2",pLAr,Pwire2[igrid][iwire],OpWireSurface2[igrid][iwire]); 
      WireSurface3 = new G4LogicalBorderSurface("WireSur3",pLAr,Pwire3[igrid][iwire],OpWireSurface3[igrid][iwire]); 
      WireSurface4 = new G4LogicalBorderSurface("WireSur4",pLAr,Pwire4[igrid][iwire],OpWireSurface4[igrid][iwire]); 
      OpWireSurface1[igrid][iwire]->SetMaterialPropertiesTable(JC_mt);
      OpWireSurface2[igrid][iwire]->SetMaterialPropertiesTable(JC_mt);
      OpWireSurface3[igrid][iwire]->SetMaterialPropertiesTable(JC_mt);
      OpWireSurface4[igrid][iwire]->SetMaterialPropertiesTable(JC_mt);
    }
  }


  return world_phys;

}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

// specific method to G4UserLimits:= SetUserMinEkine
void DMXDetectorConstruction::SetRoomEnergyCut(G4double val)
{
  // set minimum charged particle energy cut - NB: for ROOM
  theRoomMinEkine = val;
  if (theUserLimitsForRoom != 0) 
    {
      theUserLimitsForRoom->SetUserMinEkine(val); 
      G4cout << " Changing Room energy cut to: " << G4BestUnit(val,"Energy")
	     << G4endl;
    }
}  

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

// specific method to G4UserLimits:= SetUserMinEkine
void DMXDetectorConstruction::SetEnergyCut(G4double val)
{
  // set minimum charged particle energy cut - NB: for Xenon Detector
  theMinEkine = val;
  if (theUserLimitsForDetector != 0) 
    {
      theUserLimitsForDetector->SetUserMinEkine(val);
      G4cout << "Changing Detector energy cut to: " << G4BestUnit(val,"Energy")
	     << G4endl;
    }
}  

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

// specific method to G4UserLimits:= SetUserMaxTime
void DMXDetectorConstruction::SetRoomTimeCut(G4double val)
{
  // set room time cut:
  theRoomTimeCut = val;
  if (theUserLimitsForRoom != 0) 
    {
      theUserLimitsForRoom->SetUserMaxTime(val);
      G4cout << " Changing Room Time cut to: " << G4BestUnit(val,"Time")
	     << G4endl;
    }
}  

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

// specific method to G4UserLimits:= SetUserMaxTime
void DMXDetectorConstruction::SetTimeCut(G4double val)
{
  // set detector time cut:
  theMaxTimeCuts = val;
  if (theUserLimitsForDetector != 0) 
    {
      theUserLimitsForDetector->SetUserMaxTime(val);
      G4cout << " Changing Detector Time cut to: " << G4BestUnit(val,"Time")
	     << G4endl;
    }
}  

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

//void DMXDetectorConstruction::UpdateGeometry()
//{
//  G4RunManager::GetRunManager()->DefineWorldVolume(Construct());
//}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....


